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ABSTRACT 

Quasiclassical approach and geometric optics allow to describe rather accurately whispering gallery modes in 
convex axisymmetric bodies. Using this approach we obtain practical formulas for the calculation of eigenfre- 
quencies and radiative Q-factors in dielectrical spheroid and compare them with the known solutions for the 
particular cases and with numerical calculations. We show how geometrical interpretation allows expansion of 
the method on arbitrary shaped axisymmetric bodies. 
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1. INTRODUCTION 

Submillimeter size optical microspheres made of fused silica with whispering gallery modes (WGM) 1 can have 
extremely high quality-factor, up to 10 10 that makes them promising devices for applications in optoelectronics 
and experimental physics. Historically Richtmyer 2 was the first to suggest that whispering gallery modes in 
axisymmetric dielectric body should have very high quality-factor. He examined the cases of sphere and torus. 
However only recent breakthroughs in technology in several labs allowed producing not only spherical and 
not only fused silica but spheroidal, toroidal 3,4 or even arbitrary form axisymmetrical optical microcavities 
from crystalline materials preserving or even increasing high quality factor. 5 Especially interesting are devices 
manufactured of nonlinear optical crystals. Microresonators of this type can be used as high-finesse cavities for 
laser stabilization, as frequency discriminators and high-sensitive displacement sensors, as sensors of ambient 
medium and in optoelectronical RF high-stable oscillators. (See for example materials of the special LEOS 
workshop on WGM microresonators 6 ). 

The theory of WGMs in microspheres is well established and allows precise calculation of eigenmodes, radiative 
losses and field distribution both analytically and numerically. Unfortunately, the situation changes drastically 
even in the case of simplest axisymmetrical geometry, different from ideal sphere or cylinder. No closed analytical 
solution can be found in this case. Direct numerical methods like finite elements method are also inefficient 
when the size of a cavity is several orders larger than the wavelength. The theory of quasiclassical methods of 
eigenfrequencies approximation starting from pioneering paper by Keller and Rubinow have made a great progress 
lately. 8 For the practical evaluation of precision that these methods can in principal provide, we chose a practical 
problem of calculation of eigenfrequencies in dielectric spheroid and found a series over angular mode number 
I. This choice of geometry is convenient due to several reasons: 1) other shapes, for example toroids 4 may be 
approximated by equivalent spheroids; 2) the eikonal equation as well as scalar Helmholtz equation (but not the 
vector one!) is separable in spheroidal coordinates that gives additional flexibility in understanding quasiclassical 
methods and comparing them with other approximations; 3) in the limit of zero eccentricity spheroid turns to 
sphere for which exact solution and series over I up to ^~ 8 / 3 is known. 9 

The Helmholtz vector equation is unseparable 10 in spheroidal coordinates and no vector harmonics tangential 
to the surface of spheroid can be build. That is why there are no pure TE or TM modes in spheroids but only 
hybrid ones. Different methods of separation of variables (SVM) using series expansions with either spheroidal 
or spherical functions have been proposed. 11-13 Unfortunately they lead to extremely bulky infinite sets of 
equations which can be solved numerically only in simplest cases and the convergence is not proved. Exact 
characteristic equation for the eigenfrequencies in dielectric spheroid was suggested 14 without provement that if 
real could significantly ease the task of finding eigenfrequencies. However, we can not confirm this claim as this 
characteristic equation contradicts limiting cases with the known solutions i.e. ideal sphere and axisymmetrical 
oscillations in a spheroid with perfectly conducting walls. 15 



Nevertheless, in case of whispering gallery modes adjacent to equatorial plane the energy is mostly concen- 
trated in tangential or normal to the surface electric components that can be treated as quasi-TE or quasi-TM 
modes and analyzed with good approximation using scalar wave equations. 

Using quasiclassical method we deduce below the following practical approximation for the eigenfrequencies 
of whispering gallery modes in spheroid: 
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where a and b are equatorial and polar scmiaxises, k = ^- is the wavenumber, ( > 1, p = ! - \m\ = 0, 1,2,... 
and q = 1,2, 3,... are integer mode indices, a q are the q-th. roots of the equation Ai(— a q ) = (Ai(z) is the Airy 
function), n is refraction index of a spheroid and x = 1 for quasi-TE and \ = 1/n 2 for quasi-TM modes. 

2. SPHEROIDAL COORDINATE SYSTEM 

There are several equivalent ways to introduce prolate and oblate spheroidal coordinates and corresponding 
cigenfunctions. 16-18 The following widely used system of coordinates allows to analyze prolate and oblate 
geometries simultaneously: 

z=^(£ 2 - S )(l-7? 2 )] 1/2 cos(0) 
y=^[(e-s)(l-V 2 )} 1/2 sin(<P) 

z = ^V, (2) 

where we have introduced a sign variable s which is equal to 1 for the prolate geometry with £ £ [1, oo) determining 
spheroids and n e [—1,1] describing two-sheeted hyperboloids of revolution (Fig.l, right). Consequently, s = — 1 
gives oblate spheroids for £ e [0,oo) and one-sheeted hyperboloids of revolution (Fig.l, right), d/2 is the 
semidistance between focal points. We are interested in the modes inside a spheroid adjacent to its surface in 
the equatorial plane. It is convenient to designate a semiaxis in this plane as a and in the z-axis of rotational 
symmetry of the body as b. In this case d 2 /4 = s(b 2 — a 2 ) and eccentricity e = \J\ — (a/b) 2s . The scalar 
Hclmholtz differential equation 

where c = kd/2 is separable. The solution is ^ = R m i(c, £)S m i(c, "q)e lm ^ where radial and angular functions are 
determined by the following equations: 
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Here A m ; is the separation constant of the equations which should be independently determined and it is a 
function on m, / and c. With substitution £ = 2r/d the first equation transforms to the equation for the 
spherical Bessel function ji(kr) if d/2 — ► in which case the second equation immediately turns to the equation 
for the associated Legendre polynomials P^iv) with X m ,i = 1(1 + 1). That is why spheroidal functions arc 
frequently analyzed as decomposition over these spherical functions. 

The calculation of spheroidal functions and of \ m i is not a trivial task. 19 ' 20 The approximation of spheroidal 
functions and their zeros may seem more straightforward for the calculation of eigenfrequencies of spheroids, 
however we found that another approach that we develop below gives better results and may be easily generalized 
to other geometries. 



Figure 1. Graphical representation of prolate and oblate spheroidal coordinate systems (£, rj, <f>). 



3. EIKONAL APPROXIMATION IN SPHEROID 



The eikonal approximation is a powerful method for solving optical problems in inhomogcneous media where 
the scale of the variations is much larger than the wavelength. It was shown by Keller and Rubinow 7 that it 
can also be applied to eigenfrcqucncy problems and that it has very clear quasiclassical ray interpretation. It 
is important that this quasiclassical ray interpretation requiring simple calculation of the ray paths along the 
geodesic surfaces and application of phase equality (quantum) conditions gives precisely the same equations as 
the eikonal equations. Eikonal equations allow, however, to obtain more easily not only eigenfrequencies but 
field distribution also. 

In the eikonal approximation the solution of the Helmholtz scalar equation is found as a superposition of 
straight rays: 



The first order approximation for the phase function S called eikonal is determined by the following equation. 



where e is optical susceptibility. For our problem of searching for eigenfrequencies e does not depend on coor- 
dinates, e = n 2 inside the cavity and e = 1 - outside. Though the eikonal can be found as complex rays in the 
external area and stitched on the boundary as well as ray method of Keller and Rubinow 7, 8 ' 21 can be extended 
for whispering gallery modes in dielectrical bodies in a more simple way. 22 To do so we must account for an 
additional phase shift on the dielectric boundary. Fresnel amplitude coefficient of reflection 23 : 
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gives the following approximations for the phase shift for grazing angles: 
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Nevertheless, direct use of this phase shift in the equations for internal rays as suggested in 22 leads to incorrect 
results. The reason is a well known Goos-Hanchen effect - the shift of the reflected beam along the surface. The 



beams behave as if they are reflected from a fixious surface hold away from the real boundary at oy = 2 kn cos e • 
That is why we may substitute the problem for a dielectric body with the problem for an equivalent body 
enlarged on a T with the totally reflecting boundaries. The parameters of equivalent spheroid are marked below 
with overbars. 

The eikonal equation separates in spheroidal coordinates if we choose S = Si(£) + 5*2(77) + S3 (</>): 
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After immediate separation of = /J we have: 



Introducing another separation constant v we finally obtain solutions: 
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where a = r] 2 = I — fi 2 /v 2 . It is now the time to turn to the quasiclassical ray interpretation 7,8 ' 29-31 of 

whispering gallery modes. The eikonal equation describes rays spreading inside a cavity along straight lines and 
reflecting from boundaries. For the whispering gallery modes the angle of reflection is close to tv/2. The envelope 
of these rays forms a caustic surface which in case of spheroid is also a spheroid determined by a parameter 
£ c . The rays are the tangents to this internal caustic spheroid and follow along eodesic lines on it. In case of 
ideal sphere all the rays of the same family lie in the same plane. However, even a slightest eccentricity removes 
this degeneracy and inclined closed circular modes which should be more accurately called quasimodes 24 are 
turned into open-ended helices winding up on caustic spheroid precessing, 25 and filling up the whole region as 
in a clew. The upper and lower points of these trajectories determine another caustic surface with a parameter 
7] c determining two-sheeted hyperboloid for prolate or one-sheeted hyperboloid for oblate spheroid. The value 
of rjc has very simple mechanical interpretation. The rays in the eikonal approximation are equivalent to the 
trajectories of a point-like billiard ball inside the cavity. As axisymmetrical surface can not change the angular 
momentum related to the z axis, it should be conserved on the ray (L z = p 2 (f> = const, p 2 = x 2 + y 2 ) as well as the 
kinetic energy (velocity) . That is why r\ c is simply equal to the sine of the angle between the equatorial plane and 
the trajectory crossing the equator and at the same time it determines the maximum elongation of the trajectory 




Figure 2. Spheroidal optical resonator as a centrifugal billiard (after 50 reflections) 

from the equator plane. Together with the simple law of reflection on the boundaries (angle of reflectance is equal 
to the angle of incidence i.e normal component is reversed at reflection 30 r r = \ r — 2a{aYi), where a is the outward 
normal to the surface unit vector. The so-called billiard theory in 2D and 3D is extremely popular these days 
in deterministic chaos studies. This theory describes for example ray dynamics and Kolmogorov-Arnold-Moser 
(KAM) transition to chaos in 2D deformed stadium-like optical cavities and in 3D strongly deformed droplets. 31 
In this paper, however, we are interested only in stable whispering gallery modes which are close to the surface 
and equatorial plane of axisymmetric convex bodies. Axisymmetric 3D billiard is equivalent to 2D billiard in (p, 
z) coordinates. In these coordinates 3D linear rays transform into parabolas and a ball behaves as if centrifugal 
force acts on it. 30 Fig 2 shows how segments of geometric rays (turned into segments of parabolas) fill the 
volume between caustic lines in a spheroid with b/a = 0.6. If all the rays touch the caustic or boundary surface 
with phases that form stationary distribution (that means that the phase difference along any closed curve on 
them is equal to integer times 2tt), then the cigenfunction and hence eigenfrequency is found. 

To find the circular integrals of phases kS (13) we should take into account the properties of phase evolutions 
on caustic and reflective boundary. Every touching of caustic adds 7r/2 (see for example 8 ) and reflection adds 
7T. Thus for Si we have one caustic shift of 7r/2 at £ c and one reflection from the equivalent boundary surface £ s 
(at the distance a from the real surface), for S2 - two times 7r/2 due to caustic shifts at ±rj c , and we should add 
nothing for S3: 

is 
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where q = 1,2,3... - is the order of the mode, showing the number of the zero of the radial function on the 
surface, and p = I — \m\ = 0, 1, 2.... These conditions plus integrals (13) completely coincide with those obtained 
by Bykov 26 ~ 28 if we transform ellipsoidal to spheroidal coordinates, and have clear geometrical interpretation. 
The integral for Si corresponds to the difference in lengths of the two geodesic curves on r\ c between two points 



Figure 3. Segments of lines on caustic surfaces determining first quantization rule 




Figure 4. Segments of lines on caustic surfaces determining second quantization rule 



Pi = (£o Vc,<f>i) an d Pi = {£c,Vc,<f>2)- The first one goes from the caustic circle of intersection between £ c and r\ c 
along r/ c to the boundary surface £ s , reflects from it, and returns back to the same circle. The second is simply 
the arc of the circle between Pi and P2 (Fig. 3). The integral for S 2 corresponds to the length of a geodesic 
line going from Pi along £ c , lowering to — r\ c and returning to i] c at P 2 minus the length of the arc of the circle 
between Pi and P2. The third integral is simply the length of the circle of intersection of £ c and r\ c . 

These are elliptic integrals. For the whispering gallery modes when rj c <C 1 and £0 — £ c <C £ c > S2 may be 
expanded into series over r] c and ( and integrated with the substitutions of rj = r] c simp, ( = (£ 2 — £ 2 )/£ 2 - Finally, 
expressing spheroidal coordinates £ c and expressing £ through parameters of spheroid, we have: 
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Figure 5. Caustic circle determining third quantization rule 
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Now wc should solve the following system of equations: 
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Using the method of sequential iterations, starting for example from nk^a = I, Co°^ = 0, r^ ' = this system 
may be resolved: 
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where for the convenience of comparison we introduced /3 q — [§7r(g — |)] 2,/3 - The value of cos# needed for the 
calculation of $ r (10) one may estimate as cos# = yfl — (I + 1/2) 2 / '(nka) ~ y^? -1 / 3 . 
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Figure 6. Comparison of the precision of calculation of eigenfrequencies in spheroid for T\Eioo,ioo,i mode. 

The first three terms for nka were obtained in 3, 26-28 from different considerations, the last three are new. 

To test this series we calculated using finite element method (FEM) eigenfrequencies of TE modes in spheroids 
with different eccentricities with totally reflecting boundaries for I = m = 100 (Fig. 6). Significant improvement 
of our series is evident. The divergence of the series for large eccenricities is explained by the fact that the 
approximation that we used to calculate the integrals (15) but not the method itself breaks down in this case. 
Namely r/ c becomes comparable to £ and should not be treated as a small parameter. If we put a = 6, then all six 
terms in the obtained series coincide with that obtained in 9 from exact solution in sphere with a difference that 
Airy function zeros a q ~ (2.3381,4.0879,5.5206,...) stand in Schiller's series instead of approximate (3 q values 
(a g — j3 q ~ 0.017; 0.0061; 0.0033, ...). The reason is that the eikonal approximation breaks down on caustic, 
where its more accurate extension with Airy functions is finite. 7 To make our solution even better we may 
just formally use a q instead of q , hence obtaining final formula (1). The eikonal equation for the sphere may 
be solved explicitly and the expansion of the solution shows that quasiclassical approximation breaks down on 
a term (9(/ _1 ), and of the same order should be the error introduced by substitution of vector equations by a 
scalar ones. 

We may now calculate the dependence of mode separation on three indices with up to 0(l~ 2 precision: 
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It is interesting to note that when a = 2b (oblate spheroid with the eccentricity e = V0.75, the eigenfrequency 
separation in the first order of approximation between modes with the same I becomes equal to the separation 
between modes with different I and the same I — m (free spectral range). The difference appears only in the 
term proportional to 0{l~ 2 ^). This situation is close to the case that was experimentally observed in. 3 This 
new degeneracy has simple quasigeometrical interpretation - like in case of a sphere geodesic lines inclined to 
the equator plane on such spheroid are closed curves returning at the same point of the equator after the whole 
revolution, crossing, however, the equator not twice as big circles on a sphere but four times. 



4. ARBITRARY CONVEX BODY OF REVOLUTION 



To find eigenfrequencies of whispering gallery modes in arbitrary body of revolution one may use directly the 
results of the previous section by fitting the shape of the body in convex equatorial area by equivalent spheroid. 
In fact the body should be convex only in the vicinity of WG mode itself. For example a torus with a width 2Rt 
and a height tt may be approximated by a spheroid with a = Rt and b — \JRttt- Nevertheless, more rigorous 
approach may be developed. 

The first step is to find families of caustic surfaces. This is not a trivial task in general nonintegrable case. 
The following approximation may be used to find the first family of caustic surfaces 8 which are the place of 
biinvolute curves for characteristic for a chosen WGM mode geodesic lines on the surface: 

a c (s) ~ -I K 2 r i/3 (s) + 0(K4) 
cos 9(s) = K r~ 1/3 (s), (20) 

where <7 c (s) is the normal distance from a point s on the surface of the body to a caustic surface, k is a parameter 
of a family, and rt is the radius of curvature of the geodesic line (curvature of the surface in the direction of the 
beam) and cos 9 is angle of incidence of the beam in the point s. 

If we found a caustic surface from the first family, parametrized as p — g(z) then we can find also the second 
family parametrized as h(z) orthogonal to any surface of the first family with different n. 

A geodesic line for the surface is given by the following integral: 
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where p c — g(z m ) is the radius of caustic circle at maximum distance from equatorial plane. The length of 
geodesic line: 
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In analogous way for another geodesic line on a caustic surface from the other family p = h(z) we have 
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The third condition is: 
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With the substitution z = $£ c r), h(z) = f y/g - sy/l - rf in (25), and z = f^c g(z) = %y/£ 2 - s^/l - 
in (26) we again obtain expressions for spheroid obtained before (13,15). 



5. QUALITY-FACTOR 

As it was shown above, WGMs in axially symmetric bodies, parametrizes as p(z) are efficiently described by 
optical beam traveling close to surface geodesic line and suffering multiple reflections from a curved surface. If 
we want to calculate Q-factor of whispering gallery mode associated with surface reflection using quasiclassical 
approach we should take into account losses added at every reflection. The quality factor is determined by the 
following simple equation 1 : 
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where a denotes losses per unit path length. If the path is comprised of segments of polyline, then the length of 
each segment is L n = 2r^ cos 8, where is the radius of curvature of the geodesic line on the surface. Let power 
losses due to reflection in this segment is equal to T(9) and hence a n = T(8, rk)/L n . To account for total losses 
we should average a n over one coil of geodesic line and hence we finally obtain: 
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This very useful expression can be used to calculate Q-factors in arbitrary shaped whispering gallery res- 
onators associated not only with radiative but also with surface scattering and absorption. The required co- 
efficient T{8) may be deduced from the solution of model problem in a sphere, where rk = a and cos# = 
y/l - (I + 1/2) 2 / (kna) 2 ~ y /a^(l/2)- 1 / 3 arc constants: 
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Radiative losses per single reflection may be also calculated using quasiclassical arguments. 29 If a beam of light 
is internally reflected from a curved dielectric surface with radius of curvature Tk, then in the external region 
evanescent field exponentially decays as e - fer ™\/™ 2 si " 2 e ~ 1 ^ here r n is normal distance from the surface, however 
at a distance r n — (ns'm.8 — 1) tangential phase velocity reaches speed of light and the tail of evanescent field 
is radiated. These speculations allow to obtain quasiclassical estimate: 
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From the above equations one can calculate the Q-factor knowing T(8) and using the following expressions: 
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where p m = p(z m ) is a distance from the z-axis of the highest point of geodesic line, a c is the normal distance 
between surface and caustic surface. 

For radiative losses, however, only T{9) having exponential term essentially variates along the geodesic line 
with cos 9 maximal and pj~ minimal at z = for oblate geometry and vice versa for prolate one. 

Now we can calculate the radiative quality factor of WGMs in a spheroid. 
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Using substitution z = r] c cosip we obtain: 
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In conclusion. We have analyzed quasiclassical method of calculation of cigcnfrequencies and quality factors 
in dielectrica cavities and found that for spheroid they give rather precise results. 
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